function [varstd,setas_res,setas_time,obs4plot,date4plot] = ...
         get_setas4uncertainty(slon,slat,obsdate,data,depth,pos,...
         modelvar,setasfile)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% THIS FUNCTION RETURNS THE MODEL MONTHLY TIME-SERIES AND THE
% STANDARD DEVIATION OF MODEL RESULTS OVER 6 DAYS PRIOR TO THE OBSERVATION 
% DATE AT A SPECIFIED SITE. NOTE THAT THIS FUNCTION IS BASED ON THE SETAS
% MODEL MONTHLY OUTPUT.
% 
% INPUT  : - slon,slat      : site lon/lat position
%          - obsdate        : array of AWQ sampling date
%          - data           : matrix of data
%          - depth          : matrix of sampling depth          
%          - pos            : get either a surface or bottom awq
%                             observation value: options are either 's' for 
%                             surface or 'b' for bottom AWQ data 
%          - modelvar       : model variable of interest(string). 
%                             Available options are: 'temp' for temperature
%                             or 'salt' for salinity   
%          - setasfile      : monthly netcdf SETAS output file name(string)
%
% OUTPUT : - setas_res      : time-series of SETAS  monthly results
%          - setas2_time    : setas 
%          - varstd         : standard deviation over the past 6days prior
%                             to sampling date
%          - obs4plot       : observation value for the SETAS file 
%                             corresponding month 
%          - date4plot      : sampling date for the SETAS file 
%                             corresponding month
%
% Author: Benedicte Pasquer, IMOS/eMII (http://imos.org.au/)
% email: benedicte.pasquer@utas.edu.au
% May 2013   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%READ SETAS GRID AND TIME INFORMATION  FROM SETAS NETCDF MONTHLY FILE
ncid = netcdf.open(setasfile,'NC_NOWRITE');

varid = netcdf.inqVarID(ncid,'time');
time = netcdf.getVar(ncid,varid); 
varid = netcdf.inqVarID(ncid,'latitude');
lat = netcdf.getVar(ncid,varid); 
varid = netcdf.inqVarID(ncid,'longitude');
lon = netcdf.getVar(ncid,varid); 

%GET SPECIFIED VARIABLE DATA
switch modelvar
    case 'temp'
    varid = netcdf.inqVarID(ncid,'temp');
    temp = netcdf.getVar(ncid,varid,'double'); 
    modelvar = permute(temp,[2 1 3 4]);
    case 'salt'
    varid = netcdf.inqVarID(ncid,'salt');
    salt = netcdf.getVar(ncid,varid,'double'); 
    modelvar = permute(salt,[2 1 3 4]);

end
netcdf.close(ncid)

%DETERMINE GRID CELL COORDS FOR THE SITE CELL CENTRES.
%SETAS GRID 
if (~isempty(lon)),

    lon = lon' ; lat = lat' ;
    x_centre = lon ; y_centre = lat ;
    
end
% DETERMINE CELL CENTRE LOCATION
[m,n] = size(x_centre) ;
x0 = 2*x_centre(1,:) - x_centre(2,:) ;
y0 = 2*y_centre(1,:) - y_centre(2,:) ;
xf = 2*x_centre(m,:) - x_centre(m-1,:) ;
yf = 2*y_centre(m,:) - y_centre(m-1,:) ;
xs = ([x0 ; x_centre] + [x_centre ; xf])/2 ;
ys = ([y0 ; y_centre] + [y_centre ; yf])/2 ;
x0 = 2*xs(:,1) - xs(:,2) ;
y0 = 2*ys(:,1) - ys(:,2) ;
xf = 2*xs(:,n) - xs(:,n-1) ;
yf = 2*ys(:,n) - ys(:,n-1) ;
x_grid = ([x0 xs] + [xs xf])/2 ;
y_grid = ([y0 ys] + [ys yf])/2 ;

% DETERMINE  GRID CELL COORDS FOR SITE 
[isect,jsect] = xy2ij_curv (slon,slat,x_grid,y_grid,2) ;

% EXTRACT THE FULL TIME-SERIES (A MONTH) OF THE SELECTED VARIABLE AT THE CHOSEN DEPTH
% FROM GRID CELL

modelvar = (modelvar(jsect,isect,21:22,:)); 
modelvar = mean(modelvar,3); %AVERAGE OVER FIRST TWO DEPTH LEVELS
modelvar = squeeze(modelvar);
% CONVERTS TIME TO MATLAB FORMAT
time = length(modelvar);%TIME IN HOURS
%GET TIME FROM FILENAME
dash = regexp(setasfile,'-');
year = setasfile(dash-4:dash-1);
month = setasfile(dash+1:dash+2);
start_month = datenum(str2num(year),str2num(month),1);
m_length = time/24; %TIME IN DAYS
end_month = datenum(str2num(year),str2num(month),m_length);
% DATEINC: VECTOR OF SETAS DATES IN MATLAB FORMAT
dateinc = zeros(1,time);
dateinc(1) = start_month;

for tt = 2:time
    dateinc(tt) = addtodate(dateinc(tt-1),1,'hour');
end

% AWQ SITE SAMPLED ONCE PER MONTH 
% FIND INDEX OF AWQ SAMPLING DATE IN SETAS TIME VECTOR
didx = find(obsdate>start_month & obsdate<end_month);
if ~isempty(didx) 
    if  length(didx)>1
        warning('2 sampling dates in same month')
        didx = didx(1); %TAKE FIRST ONE ONLY
    end
% FIND TIMESTEPS FOR 6 DAYS PRIOR TO SAMPLING DATE
    didxsetas = find(dateinc >obsdate(didx)-5 & dateinc < obsdate(didx));
    date4plot = obsdate(didx);
    
    % GET SURFACE AND BOTTOM VALUE OF A SINGLE PARAMETER   
    obs4plot = get_valueSETASperiod(data,depth,didx,pos);   


% EXTRACT STD AT DATE OF INTEREST ONLY 
    varstd = std(modelvar(didxsetas(1):didxsetas(end)));
else
    date4plot = NaN; obs4plot = NaN;
    varstd = NaN;
end
setas_res = modelvar;
setas_time = dateinc;



